Binary black hole evolutions of approximate puncture initial data 
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Approximate solutions to the Einstein held equations are a valuable tool to investigate grav- 
itational phenomena. An important aspect of any approximation is to investigate and quantify 
its regime of validity. We present a study that evaluates the effects that approximate puncture 
initial data, based on skeleton solutions to the Einstein constraints as proposed by Faye et al. 

have on numerical evolutions. Using data analysis tools, we assess the effectiveness of these 
constraint-violating initial data and show that the matches of waveforms from skeleton data with 
the corresponding waveforms from constraint-satisfying initial data are > 0.97 when the total mass 
of the binary is > 40M© . In addition, we demonstrate that the differences between the skeleton and 
the constraint-satisfying initial data evolutions, and thus waveforms, are due to negative Hamilto- 
nian constraint violations present in the skeleton initial data located in the vicinity of the punctures. 
During the evolution, the skeleton data develops both Hamiltonian and momentum constraint vi- 
olations that decay with time, with the binary system relaxing to a constraint-satisfying solution 
with black holes of smaller mass and thus different dynamics. 

PACS numbers: 



I. INTRODUCTION 

With the developments of the past few years, numeri- 
cal relativity simulations of binary black hole (BBH) sys- 
tems from inspiral to merger are now feasible, almost 
routine. Most importantly, they are quickly becoming 
a potent tool to study highly relevant astrophysical phe- 
nomena. Approximations such as those provided by post- 
Newtonian (PN) theory have also proven to be valuable 
tools. They have the appeal of avoiding the computa- 
tional complexities associated with finding exact solu- 
tions to the Einstein field equations. As the demand 
for more efficient simulations increases, it is desirable to 
consider approximate methodologies in conjunction with 
numerical relativity approaches. A natural "marriage" 
in this regard, which is the focus of this work, is to con- 
sider full Einstein evolutions of approximately constraint- 
satisfying initial data. 

In general relativity, constructing initial data requires 
solving the Einstein constraints, a coupled set of elliptic 
equations (see Baumgarte and Shapiro 0] for a review 
on the mathematical foundations of numerical relativ- 
ity and Cook Q for constructing initial data). Thus, 
in general obtaining solutions to the Einstein constraints 
necessitates solving elliptic equations, which is a com- 
plex numerical problem. When black hole (BH) excision 
is used, the solvers are non-trivial d, [H, @ because of 
the excision boundaries. Even without excision, develop- 
ing constraint solvers is demanding Q and often requires 
introducing simplifying assumptions such as spatial con- 
formal flatness. 



Flexibility is also a very important issue. The family 
of problems addressed by numerical relativity is quickly 
expanding, involving non-traditional BH systems beyond 
the two-body problem d, Q . Without modifications to 
the standard initial data methodology, there will be lim- 
itations on the class of problems one is able to consider. 

The focus of the present work is on the full Einstein nu- 
merical evolution of constraint-violating or approximate 
initial data. Evolutions of constraint-violating BBH ini- 
tial data have been considered in the past. They were 
mostly done in the context of superposed Kerr-Schild 
BHs HH, EH, [U[ ■ More recently, constraint- violating 
initial data for punctures has been used for multiple BH 
evolutions [1,0. 

The difference with previous studies is in the building 
blocks used to construct the data. In Refs. [1, [Til ], 
the initial data sets were built from perturbative solu- 
tions of single punctures (boosted and/or spinning). Our 
approach, on the other hand, follows closely the skele- 
ton solutions of the Einstein equations introduced by 
Faye et al. These solutions are derived from the 

full Arnowitt-Dcser-Misner (ADM) Hamiltonian with the 
BHs represented by point-like sources modeled by Dirac 
delta function distributions. We consider configurations 
of non-spinning, equal-mass BBHs in quasi-circular or- 
bits and investigate how well the evolution of these ini- 
tial data is able to reproduce the corresponding results 
of constraint-satisfying initial data. We assess the ef- 
fectiveness of the skeleton initial data by computing the 
matches with waveforms from constraint-satisfying initial 
data evolutions. We find that the differences in the cvo- 
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lutions, and thus waveforms, are due to negative Hamil- 
tonian constraint violations present in the skeleton initial 
data. We observe that, during the course of the evolu- 
tion, the skeleton data develops both Hamiltonian and 
momentum constraint violations which both propagate 
away and decay over time while the binary system relaxes 
to a constraint-satisfying solution with BHs of smaller 
mass and thus different dynamics. 

In Sec. UH we derive the procedure for constructing 
skeleton puncture initial data. In Sec. IIII| we focus on 
quasi-circular configurations of equal-mass, non-spinning 
BBHs, and, using the effective potential method [||, we 
compare binding energies between skeleton and corre- 
sponding constraint-satisfying initial data. In Sec. HV[ we 
investigate the structure of the Hamiltonian constraint 
violations in the skeleton data. In Sec. [V] we present 
results of the evolutions. Sec. IVII presents an analysis 
of the nature of the constraint violations with a model 
involving a single puncture. In Sec. IVII1 we discuss the 
impact of using waveforms from skeleton evolutions on 
data analysis. Conclusions arc given in Sec. IVIIII 

The numerical simulations and results were obtained 
with the MayaKranc infrastructure as described in 

Refs. El EE IS d. 

II. SKELETON INITIAL DATA 

The traditional approach to constructing initial data 
in numerical relativity involves specifying the pair 
, Kij } , where g^ is the intrinsic 3-metric to a t = con- 
stant hypersurface E t , and denotes its extrinsic cur- 
vature. We use the index convention that Latin indices 
in the first part of the alphabet denote 4-dimensional 
spacetime indices and those from the middle denote 3- 
dimcnsional spatial indices. The pair {gij, Kij} must 
satisfy the Einstein constraint equations: 

R + K 2 - KijK ij = 16 7T p (1) 
VjK ij - V'K = Sirf. (2) 

Equations (JT|) and are respectively known as the 
Hamiltonian and momentum constraints. The operator 
Vj denotes covariant differentiation with respect to gij 
and Rij its associated Ricci tensor. Wc follow the nota- 
tion K = g^Kij and R = g^Rij. 

Although we are interested in vacuum spacetimes of 
BH systems, we have kept the matter sources p (total 
energy density) and j l (momentum density). This is so 
we are able, as in Ref. to represent the BHs as point- 
like sources modeled with Dirac delta distributions. 

The constraints Eqs. (JTJ) and © yield four equations; 
there are, thus, eight freely specifiable pieces in the data 
{gij, K^}. These free data can be used to single out 
the physical system under consideration (e.g. orbiting 
binary BHs) as well as to simplify solving the Einstein 
constraints. An elegant approach to identify the four 



pieces in {gij, Kij} that are fixed from solutions to the 
constraints was given in York [l9l |. based on work by 
Lichnerowicz [2{| and others. The method is based on 
the following conformal transformations and tensorial de- 
compositions: 

9ij = ^Qij (3) 
Kij = Aij + ^gijK (4) 

A ij = iP~ w A ij (5) 
K = K (6) 
A ij = Al j + (I,W) ij , (7) 

where A 1 i = A 1 j = and Vjj4* = with Vj covariant 
differentiation with respect to the conformal metric g^. 
In the tensorial decomposition of A % * given by Eq. {7J, 
A % i gives the transverse part of A 1 ^ , with the longitudinal 
part given by 

(LW)ii = 2 V (i W j} - 1 9i j V k W k . (8) 

With the transformations Eqs. ([3][7]), the constraint 
Eqs. (Q} and ^ become: 

8 Aip - V R - r,^ 5 K 2 + ^r 7 A %3 A l] = -16 tt^ 5 p(9) 

(A L W) i -^ip- (> V i K = 8TnJj w j i , (10) 

with R the Ricci scalar associated with the conformal 
3-metric and (A L W) 1 = Vj(hW) ij . 

At this point, we introduce the assumptions of con- 
formal flatness g^ = rjij and vanishing of both K and 
A l £ . These assumptions exhaust the eight freely specifi- 
able conditions at our disposal on {gij, K^}; five are in 
g^ , one in K and two in A % £ . The constraints then take 
the form: 

Alp + - i>~ 7 (hW) 2 = -2 mP 5 p (11) 
8 

(A L WY =8 7r^°~f , (12) 

where (LIT) 2 = (l.W) ij (l.W)ij . In the absence of mat- 
ter sources, or if one sets f = ip 10 ^, the constraints 
Eqs. (fTTj) and ([T2^1 decouple. That is, one can solve first 
Eq. (fT2|) for W l and use this solution to solve Eq. (fTTj) 
for tp. 

Following Ref. [l| , with the help of the momentum con- 
straint Eq. (fl"2"|) , we notice that 

(LIT) 2 = 2(hW) ij ViWj 

= 2 V, [(LIT) ij Wj] — 2 WjVi (LIT) ij 

= 2V l [{UV) ij Wj]-lQ-K^ w Wfi ] . (13) 

Substitution of Eq. (|13|) into the Hamiltonian constraint 
Eq. Jill) yields 

Aip + - ip- 7 Vi[(I.W) i: >Wj] = -2tt[V; 5 / 5~ ip 3 Wif] . (14) 
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We address now the matter sources. The stress-energy 
tensor for a set of non-interacting point-like particles with 
rest mass Ma, 4-velocity U A , and comoving number den- 
sity Ha is given by 



T ab = J2MaN a U%U, 



A ■ 



(15) 



where the sum is understood to run over all the particles. 
For each particle A located at x\, the comoving number 
density is given by a 5-function as 



1 



(r)]d7 



—-—S^x'-xiit)] 
aU\y/g 

5a 



iA\fg 



(16) 



with ( 4 ) g the determinant of 4-dimensional spacetime 
metric, 6 a = S 3 (x l — x A ), a the lapse function, ^ A = 
aU A = —N a U A , and N a the future-directed unit normal 
to the hypersurface St. The stress-energy tensor can then 
be rewritten as 



T 



ab 



A 



M A S A 
1a Vg 



U a A U b A . 



(17) 



Bowen and York [2l| found a solution to the momen- 
tum constraint as given by Eq. (|21[) . The solution rep- 
resent BHs with linear momentum P\ and is explicitly 
given as 



E^( 7pi 



+ n l njP j ) 



(23) 



with n l the unit normal of constant r spheres in flat 
space. In terms of Eq. |23|) . (hW) 13 takes the form: 



(hw) ij = £ ^ 

A LT 



— n l n 3 ) rik P 



In Eqs. (J23) and (ED), r A = \ \x l 



(x i 



(24) 
-x l A )/r A 



with x\ the coordinate location of BH^ . It can be shown 
that the total ADM linear momentum is P l = J2a ^a ■ 

We now turn our attention to the Hamiltonian con- 
straint Eq. (|2"D|) . As pointed out in Ref. [l[, the term 
ip 7 Vi[(l.W) ij Wj] in Eq. (JM]) is a "flesh" term that pro- 
vides the field between the particles and has the following 
contribution to the Hamiltonian: 



Given Eq. (|17p . the matter sources take the form: 
N a N b T ab 



P 



y Ma 1a 5a 

A 



(18) 



and 



j_a NcT bc 



E 



M A ±i U b A 5 A 



^P%5a 
A 



^ V 10 



(19) 



where we have used y/g = ?/> 6 y/rj, g a b 
and_L£= {i) g ac g cb . In deriving Eq. ([T9 



g ab + N a N b 



= ( 4 ) 

we have also in- 
troduced the spatial momentum vector P A = Ma^ -L b 
U A . The vector P a is related to the spatial part of the 
4- momentum p a = MU a of the point-like particles by 
P a = ip 4 ±f p b . Substitution of the source Eqs. JT8J) 
and JTSJ) into Eqs. O and fTQ yields 



Ai/j + -XfVi^Wf^Wj] = -2tt 



Vv 



(20) 



where 



m A 



Mai A WiP A 



^ 7 



(21) 



(22) 



^i[(hW) ij Wj]d 3 x = -7 J ^{UV) l3 W J V t ipd i x . 



The only approximation that goes into defining the skele- 
ton initial data is to neglect the contribution from this 
term. With this approximation, the Hamiltonian con- 
straint Eq. (|2"D|) reads: 



m A S A 



(25) 



with m A given by Eq. (|22p . Notice that m A is singular 
at x % = x\ because ip and W % are singular at x A . Follow- 
ing Ref. [l|, we solve Eq. l[^5l) by means of Hadamard's 
"partie finie" procedure [22j |: that is, 







1-4^A- X 
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m A {x l )5 A \ 
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(26) 
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where 



,(rcg) _ 



Mai a WfPj 

<f>A 



&A 



111 



(reg) 



1A 



W A P\ 



1 



B^A 
P l P, 



2tab 

1/2 



v / -1 



[IPs Pi 



(n? B P A )(n AB Ph)\ 



(27) 
(28) 

(29) 
(30) 



with tab 



and 



L AB 



,)/r A B- The 



parameter vtla is commonly known as the bare mass of 
the BH. On the other hand. M. is known as the irre- 
ducible mass of the BH. $a is the regularized value of 
tp(x A ). In summary, the skeleton initial data {gij, Kij} 
is then given by g^ = ip 4 Vij anc i Kij = 2 (ILV^)^- with 



tp given by Eq. ([25)1 and (LW) 4J given by Eq. 

For comparison, the exact or constraint-satisfying 
puncture initial data method [23| consists also of = 
ipSij and Kij = %l)- 2 (LW)ij with (LW) ij given by 
Eq. ([24]) . but in this case 



A 2TA 



with u a regular solution to 



Au- 



1 



1 ip 



0. 



(31) 



(32) 



III. QUASI-CIRCULAR INITIAL DATA 



We restrict our attention to initial data configurations 
representing two equal mass (A4i = M2 = -M, toJ 108 ' = 
m (reg) _ 772) , non-spinning BHs in quasi-circular orbits. 



That is P[ = -Pi = P\ r i2 = \\x 
n\ 2 P l = 0. Under these assumptions: 



where 



m m 

ib = 1 H 1 

V 2n 2r 2 



7W 7 7 P 2 



= d, and 



(33) 



4> 



4rf$ 7 



$ = 



7 



m 

1+ 2d 



P 2 



M 2 <Z> 4 



1/2 



(34) 
(35) 

(36) 



While deriving Eq. ([34]) . we used that for circular orbits 
W i P i = TP 2 /{Ad) with P 2 = piPt = /"/' , //, / as can be 
seen from Eq. (j5D|) . 



We focus now on the differences between the 
constraint-satisfying and skeleton initial data for 
quasi-circular sequences using the effective potential 
method The general idea of this method is to find 
configurations that satisfy the condition: 



dE b 
dl 



= 0, 



(37) 



M.J 



with Ef, = E — M the binding energy of the system. 
The distance I is a measure of the proper separation be- 
tween the BHs (e.g. horizon to horizon), and M = 2 Ai 
is the sum of the irreducible masses. The quantities E 
and J are respectively the total ADM mass and angular 
momentum of the system p4| . which can be computed 
from: 



E = ~ (f> V^d 2 S l 

1 7T 



Ji 



3£ (l> x^K kl d 2 Si 
8tt 



It is not too difficult to show from Eq. 



(38) 
(39) 
that, given 

= tp~ 2 (hW)ij the ADM angular momentum for bi- 
naries initially in quasi-circular orbits is J — d P. On the 
other hand, with ip given by Eq. (|3"3"|) the total ADM mass 
from Eq. (piS)) is given by the sum of the bare masses of 
the BHs, namely E = 2 m; thus, the binding energy be- 
comes Ei, = 2 m — 2 A4 . The bare masses for the skeleton 
initial data are obtained by solving the implicit Eq. (|3~4"]) 
using a Newton-Raphson method. 

Figure[T](top panel) shows the comparison of the bind- 
ing energy E^ as a function of the total ADM angular mo- 
mentum J between the constraint-satisfying initial data 
from Tichy and Brugmann [25| (squares) and the skele- 
ton initial data in this work (triangles). The lower panel 
in Fig. [1] shows the corresponding % relative difference 
between both results. Not surprisingly, as the binary 
separation increases (i.e. larger angular momentum), the 
differences diminish. For reference, the vertical lines de- 
note the angular momentum for typical data sets consid- 
ered in the literature: QCO in Ref. [26|. Rl in Ref. [27} 
and D10 in Ref. [2j|. The differences in binding energy 
between the skeleton and the constraint-satisfying initial 
data are - 20 % for QCO, - 6 % for Rl and - 2 % for 
D10. 

TableUprovides the parameters of the initial configura- 
tions for both the skeleton and constraint-satisfying data 
sets. The cases of exact or constraint-satisfying initial 
data are labeled with the letter "e" and the correspond- 
ing skeleton or approximate case with the letter "a" . 

As mentioned before, the only fundamental difference 
between the two initial data sets is in the conformal factor 
-0. For the constraint-satisfying data set ip is computed 
from Eq. (|3ip by solving the Hamiltonian constraint in 
the form given by Eq. ([32)1 and for the skeleton the con- 
formal factor ip is constructed using Eq. (|2"6"1) . In Fig. [21 
we show the relative difference Sip/tp = (ip a — ip e )/ijj e 
from the two data along the axis joining the punctures 
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FIG. 1: Comparison of the binding energy Eb as a function of 
the total ADM angular momentum J between the initial data 
from Tichy and Briigmann [25| (squares) and the skeleton 
initial data (triangles) 
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FIG. 2: The relative difference in the conformal factor 
i/j between the skeleton initial data and the corresponding 
constraint-satisfying data along the rr-axis joining the punc- 
tures for the three cases labeled in Fig. [T] 



Run 


d/M 


P/M 


m/M 


M/M 


E/AI 


J/M 2 


QCOc 


2.337 


0.33320 


0.45300 


0.519071 


1.0195 


0.7787 


QCOa 


2.337 


0.33320 


0.48950 


0.519071 


0.9790 


0.7787 


Rlc 


6.514 


0.13300 


0.48300 


0.505085 


0.9957 


0.8664 


Rla 


6.514 


0.13300 


0.49717 


0.505085 


0.9943 


0.8664 


DlOe 


10.00 


0.09543 


0.48595 


0.500000 


0.9895 


0.9530 


DlOa 


10.00 


0.09543 


0.49458 


0.500000 


0.9891 


0.9530 



TABLE I: Initial data parameters: The punctures have bare 
masses m, linear momenta ±P and are separated by a dis- 
tance d. The irreducible mass of each BH from m (reg) is M. 
The ADM masses and angular momenta of the spacetimes are 
given respectively by E and J 



( x-axis) for the D10, Rl and QCO cases. Notice the large 
differences in the immediate vicinity of the punctures. In 
the next section, we will investigate how these differences 
translate into constraint violations. 



IV. HAMILTONIAN CONSTRAINT 
VIOLATIONS 

For the remainder of the paper we will concentrate 
our attention on the D10 case: a situation in which the 
BHs are not too close to the merger and with an initial 
separation that permits a reasonable overlap with the 
post-Newtonian regime (28l . [2^ | . It is important to point 
out that the numerical data DlOe, although called ex- 
act, also violate the constraints initially. The violations 



in the exact initial data, however, are a consequence of 
numerical errors which can be made arbitrarily small in 
the limit to the continuum. On the other hand, the con- 
straint violations in the skeleton data arc independent of 
the resolution of the computational grid. 

In order to understand the nature of the constraint vi- 
olations in the skeleton initial data and in particular their 
dynamics in the course of the evolution, we take the point 
of view that the violations introduce "spurious" sources 
p and j l in Eq. (p} and @ , respectively. Notice that ini- 
tially we do not have a "spurious" momentum density j 1 
because the skeleton initial data by construction are an 
exact solution to the momentum constraint. It is impor- 
tant to keep in mind that one should not assign physical 
properties to p and j l . They are only used to quantify 
constraint violations. In particular, the violations p are 
not restricted to satisfy energy conditions and thus are 
free to take negative values. 

Fig. [3] shows a surface plot of p for the BBH skele- 
ton initial data in the neighborhood of one of the punc- 
tures. Notice that the puncture seems to be embedded 
in a "cloud" or a pocket of negative p. Furthermore, 
the cloud is more negative in the direction aligned with 
the linear momentum of the puncture (in this case the 
y-axis). This effect is more evident from Fig. 2] where 
we plot p in the top panel along the x-axis (the direc- 
tion joining the BHs) and in the bottom panel along the 
y-axis. The glitches at the bottom of the bottom of the 
constraint violation pockets are due to refinement bound- 
aries. 
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FIG. 4: Sources p corresponding to Fig. [3] along the s-axis 
joining the BHs (top panel) and along the j/-axis (bottom 
panel), the linear momentum direction. 



V. SKELETON EVOLUTIONS 

Given the initial data, we turn our attention now to 
evolutions. The evolution runs were done on a compu- 
tational grid with 9 refinement levels, the finest 5 levels 
containing 24 3 gridpoints in radius and the remaining 4 
with 48 3 gridpoints in radius. To check the dependence 
of the results with resolution, we considered grid spac- 
ings at the finest level of M/38.4, M/44.8 and M/51.2. 
The results presented here were done at the resolution of 
M/51.2. 

Fig. [5] shows the trajectory of one of the BHs from the 
skeleton initial data (dashed line) as well as its constraint- 



FIG. 5: Trajectory of one of the BHs from the skeleton initial 
data (dashed line) as well as its constraint-satisfying counter- 
part (solid line). 

satisfying counterpart (solid line). Both trajectories are 
very close to each other during the first quarter orbit. 
Beyond that point, the BH from the skeleton initial data 
follows an eccentric orbit. Finally, near merger or at 
the plunge, the trajectories once again lie very closely 
together. 

In Fig. [6l we compare the waveforms of the skeleton 
initial data with its constraint-satisfying counterpart as 
detected at 50 M. The presence of a phase shift between 
the two waveforms is evident. The constraint-satisfying 
initial data evolution reaches the merger approximately 
10 M before the skeleton initial data evolution. This 
difference remains within 1 M of this between different 
resolutions. Another difference in the two evolutions is 
in the inspiral. As mentioned before, the skeleton data 
yields a larger eccentricity in the inspiral. This can be 
clearly observed from Fig. [JJ where the same compari- 
son as in Fig. [6] is shown but in terms of the amplitude 
(top panel) and phase (bottom panel). Here we have ap- 
plied a time shift of 10 M to align the point at which 
the waveforms reach their maximum values. The inspi- 
ral and plunge of the binary is before the "knee" in the 
phase or the maximum in the amplitude. On the other 
hand the quasi-normal ringing of the final BH takes place 
after the knee in the phase and the maximum in the am- 
plitude. Notice that the phases are practically identical 
for both cases. Furthermore, both the post-knee phase 
and post-maximum amplitude are almost the same for 
skeleton and constraint-satisfying evolutions, which is an 
indication that the final BHs are almost identical [3(| ■ On 
the other hand, the inspiral amplitudes in Fig. [JJ clearly 
show differences in the level of eccentricity as seen by the 
oscillations in the amplitude. 
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FIG. 6: Real parts of the waveform, r ^ 2 ^ 2 M, extracted 
at r a = 50 M for both the skeleton (dashed line) and the 
constraint-satisfying (solid line) initial data. 



FIG. 7: Amplitude (top panel) and phase (bottom panel) 
of the waveforms T$l\ M in Fig. [6] skeleton data (dashed 
line) and constraint-satisfying data (solid line) . The time axis 
has been shifted by 10 M to align the point at which the 
amplitudes reach their maximum values. 



From the waveforms, we have computed the energy 
"^rad anc ^ angular momentum </ rad radiated. For the 
constraint-satisfying initial data, we obtained -E rad = 
0.0354 ill and J rad = 0.3060 M 2 and for the skeleton 
data £ rad = 0.0359 M and J rad = 0.3063 M 2 , which 
correspond to differences of 1.4% and 0.1% respectively. 
These differences arc consistent with differences in energy 
and angular momentum in the initial data (< 10 _4 %). 

To better understand the change in trajectories and 
the corresponding phase shift reflected in the waveforms 
(see Fig. [5]), we have tracked the evolution of the ap- 
parent horizon (AH) masses. The AH mass for one of 
the BHs is plotted in Fig. [5] where the error due to grid 
spacing resolution is of order 10~5M. While the AH 
mass for the constraint-satisfying evolution stays rela- 
tively constant (solid line) , the AH mass for the skeleton 
evolution varies significantly (dashed line). In fact, the 
mass starts more than 10% higher than the constraint- 
satisfying value and monotonically decreases. Empiri- 
cally, the AH masses decrease as 1/t at late times. By 
fitting a polynomial in 1/t to the AH evolution at late 
times, we find the mass asymptotes to 0.501 ± 0.001M, 
within 0.2% of the constraint-satisfying initial AH mass. 
However, the BHs merge before the skeleton AH mass 
could reach this asymptotic value. The differences and 
evolution of the AH masses early in the evolution of the 
skeleton data are consistent with the picture of a binary 
whose masses and therefore binding energy and dynamics 
are altered. 



VI. SINGLE PUNCTURE ANALYSIS 

As noted in Sec. IIII1 the Hamiltonian constraint viola- 
tions are negative in the vicinity of the punctures. To bet- 
ter understand the evolutions of the skeleton initial data, 
we consider a test case where we evolve a single, non- 
spinning puncture and add by hand negative constraint 
violations surrounding it. That is, we solve the Hamilto- 
nian constraint as if there were an additional matter field 
p present, namely 

Atp = —2n pip 5 . (40) 

In order to guarantee the existence of a solution as dis- 
cussed in fl9i . [3l| . if p > 0, one needs to re-scale the 
source according to the conformal rescaling p = pip~ s , 
with s > 5 . In our case, however, we are mostly inter- 
ested in p < 0, which does not require any rescaling for 
existence of a solution. 

Following the procedure for multiple BHs, see Eq. (|3ip . 
we use the ansatz ip = ip + it, with ip a = 1 + m/2r the 
solution to the homogeneous equation (i.e. the single 
puncture solution). With this ansatz, Eq. (|4TJ|) becomes 

Au = -2irp(ijj + u) n (41) 

where n = —3 for p > and n = 5 for p < 0. 
For simplicity, we choose 

p = ^iTe-O-^) 2 /^ 2 ) (42) 

where ro is the position with respect to the puncture and 
rn = for p > and m = —5 for p < 0. The factor ip™ is 
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FIG. 8: The evolution of the AH mass of one of the BHs shown 
for both the constraint-satisfying initial data DlOe (solid line) 
and its skeleton counterpart (dashed line). Errors in AH due 
to grid spacing are of order 10 5 M. 



Model 


F/M 2 


E/M 


M\ H /M 


M f AH /M 


E p /M 


Fi 


0.001 


1.0046 


1.0012 


1.0041 


0.0034 


F 2 


-0.001 


0.9902 


0.9973 


0.9911 


-0.0071 


F 3 


-0.010 


0.9102 


0.9858 


0.9183 


-0.0756 



TABLE II: Models: Results of evolutions a single puncture 
in the presence of a Gaussian source p with r = a = 1 M 
and amplitude F . The initial AH mass and ADM energy are 
M\n and E respectively. The AH mass at the end of the 
simulation is given by M S AH /M, and E P = E- M l AH . 



necessary for regularity of the solution u at the puncture. 
We also assume that the source p does not have initial 
momentum (i.e. f = 0); thus, the momentum constraint 
remains satisfied as in the vacuum case. 

Table [TT] lists the results from the evolutions for r D = 
a = 1 M. Notice that model F± has a positive source (i.e. 
F > 0) and the other two have negative Hamiltonian con- 
straint violations. The effect of the source p is evident 
in the ADM mass (E) and initial AH mass (M l AH ). For 
the positive source, the masses are larger than the punc- 
ture mass in vacuum, 1 M, and smaller for the negative 
sources. Also in Table HH we include E p = E — M. AH , 
which gives a measure of the extra energy content in the 
initial data due to p. 

We evolved the models in Tabic HQ for 300 M. Fig. E 
shows how the AH mass evolves during the evolution. 
We have evolved the model F3 at different resolutions 
and estimated the AH masses to have an approximate 
relative error of ~ 0.002%. We observed that at late 



FIG. 9: Evolution of the AH mass for the models described 
in Table HIl The error between resolutions for F3 was of order 
10~ 5 over the course of the evolution. 

times the AH mass evolves as M. AH + C/t. The values 
reported in Tabic HIl are those extrapolated to t — * 00. 

The evolutions of the single puncture models clearly 
demonstrate that depending on the signature of p, the 
mass of the BH, as measured from the AH, will increase 
or decrease. That is, over the course of the evolution, the 
AH masses evolve to approach the ADM energy, decreas- 
ing for a negative p and increasing for positive p. In other 
words, the source p initially hovering near a puncture will 
fall into the BH, increasing or decreasing its mass as the 
system becomes stationary depending on the sign of p. 
The extent to which the final mass of the BH approaches 
the total ADM energy depends on how much of the den- 
sity p is "accreted" by the BH. Since in our case, we do 
not impose the restriction of positivity on p, the BH is 
free to decrease its mass. Notice also that the final AH 
mass does not satisfy the condition M AH = M\h + E p , 
which means that a fraction (< 1% in our cases) of E p 
mass is radiated away. The choice of centering the Gaus- 
sian at r = 1 M was aimed at favoring the amount p 
accreted by the BH. 

Figure [10] shows the Hamiltonian constraint violation 
p near the beginning of the simulation at t = 0.078 M 
(top panel) and at the end, t = 300 M, of the simula- 
tion (bottom panel). Solid lines represent the constraint- 
satisfying case and dashed lines the F3 model. Fig. [IT] 
shows the corresponding results for the momentum con- 
straint violations j 1 along the x-axis. By construction, 
initially there are only Hamiltonian constraint violations 
in the F3 model. However, it is evident from the top 
panel in Fig. [11] that constraint violations in the mo- 
mentum constraint develop also very early in the evo- 
lution. The growth of momentum constraint violations 



a 

x/M 



2 
x M 



FIG. 10: Hamiltonian constraint violation p near the begin- 
ning of the simulation at t — 0.078 M (top panel) and at the 
end. t = 300 M, of the simulation (bottom panel). Solid lines 
represents the constraint-satisfying case and dashed lines the 
F3 model. The constraint violations still present at late times 
are due to discretization around the punctures. 



proceed up to a time t ~ 3 M. The subsequent dynamics 
of the constraint violations consists of ingoing and out- 
going waves. Because of the proximity to the puncture, 
the outgoing waves are a little bit weaker, with most of 
the constraint violations "accreted" by the BH. After ap- 
proximately t ~ 50 M of evolution, the F$ model relaxes 
to the configuration of the constraint-satisfying puncture 
and remains there as seen in the bottom panels in Figs. 1101 
and llll The final constraint violations in the system arise 
from numerical errors. 

An important aspect to point out is that although the 
constraint-violating cases relax to a constraint-satisfying 
solution, the solutions that they relax to are not necessar- 
ily the same solutions as a single puncture in a vacuum 
spacetime. The new solution satisfies the Einstein equa- 
tions but for a single puncture spacetime with a smaller 
mass. A similar situation occurs in the binary case; the 
system relaxes to a binary solution, but this solution is 
different than the vacuum case. The reason for this be- 
havior is not currently understood. 



VII. IMPACT ON DATA ANALYSIS 

Finally, we want to address the extent to which the 
waveforms from evolutions of skeleton initial data may 
be of use in exploring gravitational wave astronomy. We 
will focus on computing the matches between the skele- 
ton and the constraint-satisfying waveforms. In princi- 
ple, the match would be between the detector output, hi 



FIG. 11: Same as in Fig.fTOlbut for the momentum constraint 
violation j x . The constraint violations still present at late 
times are due to discretization around the puncture. 



and the template, h 2 . Here hi is the waveform from the 
constraint-satisfying evolution and h 2 from the skeleton 
initial data evolution. Specifically, we will c omp are the 
waveforms using the minimax match given by |32l.[33l.[3^|. 



Match 



max mm max ■ 

*0 *2 *1 



(hi\h 2 



(43) 



^{hi\hi){h 2 \h 2 ) 
where the inner product of two templates is defined by 



{hi\h 2 



4 Re 



tn 



hi(f)h* 2 (f) 
Sh(f) 



df- 



(44) 



The match is maximized over the time of arrival of 
the signal, to, and minimized/maximized over the initial 
phase, $1 and $2, of the orbit when the signal/template 
enters the LIGO band. The variable Sh(f) denotes the 
noise spectrum for which we use the initial LIGO noise 
curve [35j. The domain [/ m i n , /max] is determined by the 
detector bandwidth and the masses of our signal - set 
such that the initial frequency of the numerical waveform 
just enters the LIGO band. We have chosen to study 
the match for values of the total mass of the BBH sys- 
tem greater than 20M@ because of the limited number 
of cycles that our waveforms include. A more detailed 
description of our calculation of the minimax match is 
given in fvf\ . 

The match between the constraint-satisfying and skele- 
ton data versus mass is plotted in Fig. [T2J As the to- 
tal mass increases, the match between the waveforms in- 
creases, becoming > 0.99 at 60M Q . At such large total 
mass, the signal is dominated by the plunge and ring- 
down. Comparisons of the plunge and ring-down show 
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FIG. 12: Matches 

(see Fig. [7]) that the difference between the skeleton and 
constraint-satisfying evolution are very small. At masses 
lower than about 40M Q , the match drops below 0.97 due 
to the difference in the binary dynamics prior to merger. 
We note that our calculation of match did not maximize 
over the mass of the two waveforms. Maximizing over 
the mass would have diminished the differences between 
the two waveforms. 



VIII. CONCLUSIONS 

We have carried out a study of the evolution of skele- 
ton, puncture BBH initial data as proposed by Faye et al. 
[l[ . We focused on non-spinning punctures at initial sep- 
arations of 10 M, where the difference in binding energy 
with the constraint-satisfying initial data is < 2%. We 
showed that during the inspiral the skeleton data yields 
different dynamics; however, this difference significantly 
diminishes as the binary enters the plunge, merger and 
ring-down. 

We tested the match between the constraint-satisfying 
and skeleton data for a series of total masses between 
20M Q and 130M Q . Our results indicate that gravita- 
tional wave data analysis would have some tolerance for 
constraint-violating data, especially for those binaries in 
which the signal is plunge-merger dominated, as is the 
case of high mass BHs. We conclude that although the 
two systems were different, with one clearly violating the 
Einstein equations, the differences were not enough to im- 



pact the match statistics for the mass ranges we included 
in our study and for the number of cycles present in 
our numerical waveforms. If these systems were evolved 
starting with a larger initial black-hole separation, the 
constraint violations would be smaller; and, therefore, 
the waveforms generated could be useful for detection 
over the complete BBH mass range for initial LIGO. If, 
however, larger constraint violations are present in the 
data that drive the early BH mass' lower, the differences 
may lead to errors in parameter estimation. 

We also analyzed the impact of the Hamiltonian con- 
straint violations. We showed that the main feature of 
the skeleton data is two packets of negative constraint 
violations in front of and behind the BH, along the di- 
rection of its momentum. We conjectured that these neg- 
ative constraint violations acted as a source density that 
gets absorbed by the BHs during evolution. To test our 
conjecture, we considered a model consisting of a single, 
non-rotating puncture in which we artificially added a 
stationary Gaussian shell source that mimics the Hamil- 
tonian constraint violations in the skeleton data. The 
evolutions of this single puncture model reproduce the 
decrease in the mass of the BH observed in the evolution 
of the skeleton data. 

One remarkable aspect of our study is the ability of the 
BSSN equations and moving puncture gauges to evolve 
stably data away from the constraint surface. What is 
even more remarkable is how the evolution brings the 
data back to the Einstein constraint surface. We are 
currently investigating a broader class of solutions with 
this property. 

In summary, our numerical evolutions show that the 
skeleton initial data proposed by Faye et al. [l[ embeds 
the BHs in a "cloud" of negative constraint violations. 
These constraint violations act as a source field that when 
accreted by the BHs decreases their masses. The change 
in the masses modifies the binding energy of the binary 
and thus affects its orbital dynamics (e.g. adding ec- 
centricity) but had little affect on the match of the two 
waveforms for initial LIGO for high mass black holes. 
The observed effects will decrease as the initial binary 
separation increases. 
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